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Abstract 

Quantum-dynamical full-dimensional (15D) calculations are reported for the protonated water 
dimer (HsO^) using the multiconfiguration time-dependent Hartree (MCTDH) method. The dy- 
namics is described by curvilinear coordinates. The expression of the kinetic energy operator in this 
set of coordinates is given and its derivation, following the polyspherical method, is discussed. The 
PES employed is that of Huang et al. [JCP, 122, 044308, (2005)]. A scheme for the representation 
of the potential energy surface (PES) is discussed which is based on a high dimensional model 
representation scheme (cut-HDMR), but modified to take advantage of the mode-combination 
representation of the vibrational wavefunction used in MCTDH. The convergence of the PES ex- 
pansion used is quantified and evidence is provided that it correctly reproduces the reference PES 
at least for the range of energies of interest. The reported zero point energy of the system is con- 
verged with respect to the MCTDH expansion and in excellent agreement (16.7 cm~^ below) with 
the diffusion Monte Carlo result on the PES of Huang et al. The highly fluxional nature of the 
cation is accounted for through use of curvilinear coordinates. The system is found to inter convert 
between equivalent minima through wagging and internal rotation motions already when in the 
ground vibrational-state, i.e., T=0. It is shown that a converged quantum-dynamical description 
of such a flexible, multi-minima system is possible. 
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I. INTRODUCTION 



The protonated water-dimer HsO^ (Zundel cation) is the smallest system in which an 
excess proton is shared between two water molecules. Much effort has been devoted to this 
system due to the importance that the solvated proton has in several areas of chemistry and 
biology. In recent years a fast development of the spectroscopical techniques available to 
probe the vibrational dynamics of ionic species in the gasphase has taken place, and several 
studies have appeared around the HsO^ system [l,[2, 3, J| and other more complex clusters 
and molecules ^, 6]. In order to achieve a satisfactory understanding of the spectroscopy 
and dynamics accurate theoretical results are needed in parallel and several works have ap- 
peared to fill this gap which are based on different theoretical approaches, from classical- to 



quantum-dynamical methods 



Accurate quantum-dynamical simulations 



of the dynamics and spectroscopy of the system require of an accurate reference potential 
energy surface (PES). Very accurate PES and dipole moment surfaces (DMS) have been 
recently produced by Huang et al. 7|] for HsO^. They are based on several ten-thousands 
of coupled-cluster calculations combined with a clever fitting algorithm which is based on a 
redundant set of coordinates, namely all atom-atom distances. The PES is able to describe 
the floppy motions occurring at typical energies of excitation in the linear IR regime and 



dissociates correctly 0]. Several works have already appeared in which the PES and DMS of 
Huang et al. [3] were used 4,[ll|, [l^. In Ref. [l^ full- dimensional vibrational calculations 
for HsO^ were undertaken using diffusion Monte Carlo (DMC) and variational wavefunction 
methods. The zero point energy (ZPE), some vibrational excited-state energies and proper- 
ties of the ground vibrational state of the system are reported and discussed. The PES of 
Huang et al. was also used in Ref. ^ where the experimental IR spectrum of the system 
was reported together with the calculation of relevant excited-vibrational states. 

As will be later discussed, HsO^ features several symmetry-equivalent minimum energy 
structures and large amplitude motions between different regions connected by low energy 
barriers. For such systems curvilinear coordinates, e.g. bond-bond angle, dihedral angle, or 
internal rotation coordinates, offer an advantage over rectilinear coordinates since they pro- 
vide a physically meaningful description of the different large-amplitude molecular motions. 
As a remark, an attempt was made to describe the system by a set of rectilinear coordinates, 
since they result in a simple expression of the kinetic energy operator (KEO). The resulting 



3 



Hamiltonian was unsuccessful at describing of variuos large-amplitude displacements, which 
resulted in abandonement of such an approach and introduction of a curvilinear-coordinates 
based Hamiltonian. When deriving the KEO for a Hamitonian in curvilinear coordinates 
one may define those in terms of rectilinear atom motions and use standard differential cal- 
cuius to obtain the desired expression |13!]. This approach becomes extremely complex and 
error prone already for systems of a few atoms. Another possibility is to use an approximate 
KEO where some cross-terms are neglected, thus simplifying the derivation. In this case, 
however, an uncontrolled source of error is introduced in the calculation. These drawbacks 
are overcome by employing the polyspherical approach [isl, 16] when defining the co- 
ordinate set and deriving the corresponding KEO. In the polyspherical method the system 
is described by a set of vectors of any kind (e.g. Jacobi, valence, . . . ). The kinetic energy 
is given in terms of the variation in length and orientation of these vectors, the latter is 
defined by angles with respect to a body fixed frame. Following the polyspherical method 
the KEO of the system at hand is derived in a systematic way without use of differential 
calculus. The polyspherical approach has already been discussed in a general framework as 
well as for some molecular systems (see for instance , lid . IITI]). In this paper the full 

derivation of the KEO in a set of polyspherical coordinates is illustrated for the HsO^ cation 
and all the key steps are carefully discussed. The derivation is followed step-by-step for this 
rather complex system, thus providing the basic guidelines to be followed for the derivation 
of KEOs for complex molecular systems and clusters. 

In quantum-dynamical wavefunction-based studies, as the one being introduced here, it 
is convenient to represent the operators on a discrete variable representation (DVR) grid, 
which in the multidimensional case is the direct product of one-dimensional DVRs. The 
potential operator is then given by the value of the PES at each point of the grid. In the 
present case, however, the resulting primitive grid has more than 10^^ points. This number 
makes clear that the potential must be represented in a more compact form to make the 
calculations feasible. Algorithms exist which take advantage of the fact that correlation 
usually involves the concerted evolution of only a small number of coordinates as compared 
to the total number. Hierarchical representations of the PES are then constructed including 
potential-function terms up to a given number of coordinates fisl . fiol . 20|. We discuss here 
a variation of a hierarchical representation of the PES in which the coordinates are grouped 
together into modes, and the modes are the basic units that define the hierarchical expansion. 
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It will be shown that this approach allows for a fast convergence of the PES representation 
if the modes are defined as groups of the most strongly correlated coordinates. 

The present work provides new reference results on the properties of the ground vibra- 
tional state of HsQt. The convergence of the given results is established by comparison to 



the DMC resuhs 



12| on the reference PES of Huang et al. Q]- These results, together with 



the Hamiltonian defined here, set also the theoretical and methodological framework used 



2l| where the IR spectrum and vibrational dynamics of HsO^ are 



in the companion paper 
analyzed. 

The paper is organized as follows. Section [III presents a brief description of the multicon- 
figuration time-dependent Hartree (MCTDH) method. Section IIIII details the construction 
of the Hamiltonian operator for the HsO^ system. The derivation of the KEO is discussed 
in nil A| while the construction of the potential is detailed in IIIIBI Section IIV Al discusses 
the calculation of the ground vibrational state of the system and gives a comparison to 
available results. In section IIVBI some properties of the system in relation to its ground 
vibrational-state wavefunction are discussed. In section IIV CI the quality of the potential 
expansion is analyzed and discussed. Section |V] provides some general conclusions. 



II. BRIEF DESCRIPTION OF THE MCTDH METHOD 



The Multiconfiguration time- dependent Hartree (MCTDH) method 22|, |23|, l24, l25|] is a 
general algorithm to solve the time- dependent Schrodinger equation. The MCTDH wave 
function is expanded in a sum of products of so-called single-particle functions (SPFs). 
The SPFs, (f{Q,t), may be one- or multi-dimensional functions and, in the latter case, the 
coordinate Q is a collective one, Q = {qk, ■ ■ ■ ,qi). As the SPFs are time-dependent, they 
follow the wave packet and often a rather small number of SPFs suffices for convergence. 

The ansatz for the MCTDH wavefunction reads 

ni rip p 

jl jp K = l 

J 

where / denotes the number of degrees of freedom and p the number of MCTDH particles, 
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also called combined modes. There are SPFs for the /t'th particle. The Aj = Aj^,,,j^ 
denote the MCTDH expansion coefficients and the configurations, or Hartree-products, $j 
are products of SPFs, implicitly defined by Eq. ([T]). The SPFs are finally represented by 
linear combinations of time-independent primitive basis functions or DVR grids. 

The MCTDH equations of motion are derived by applying the Dirac-Frenkel variational 
principle to the ansatz Eq. ([1]). After some algebra, one obtains 

lAj = Y,{'^j\H\'^l)Al, (2) 

L 

^(^W = (l-p(-)) (pW)~\H)("V^"^ (3) 

where a vector notation, (/j*^''^ = ((/?['^\ ■ ■ ■ , (yji']^'')"^, is used. Details on the derivation, the 
definitions of the mean-field the density p^'^^ and the projector P^'^\ as well as more 



general results, can be found in Refs. 



23 



24 



25| . Here and in the following (except for 



Section [HT]) we use a unit system with h= 1. 

The MCTDH equations conserve the norm and, for time-independent Hamiltonians, the 
total energy. MCTDH simplifies to Time-Dependent Hartree when setting all = 1. 
Increasing the recovers more and more correlation, until finally, when equals the 
number of primitive functions, the standard method (i. e. propagating the wave packet on 
the primitive basis) is used. It is important to note, that MCTDH uses variationally optimal 
SPFs, because this ensures early convergence. 

The solution of the equations of motion requires to build the mean-fields at every time- 
step. A fast evaluation of the mean-fields is hence essential. To this end the Hamiltonian is 
written as a sum of products of mono-particle operators: 

H = j^c.\[hi^\ (4) 

(k) 

where hr operates on the K-th particle only and where the are numbers. In this case the 
matrix-elements of the Hamiltonian can be expressed by a sum of products of mono-mode 
integrals, the evaluation of which is fast. 

s p 

{<!>j\H\<!>L) = J2crl[(^M''^\y^J, (5) 

r=l K=l 

Similar expressions apply to the matrix of mean-fields (H)'^^^. For further information on 



the MCTDH method, see http://www.pci.uni-heidelberg.de/tc/usr/mctdh/ 



III. SYSTEM HAMILTONIAN 



A. Kinetic energy operator for HsO^ 

The derivation of the exact kinetic energy operator for H^O^ is based on a polyspherical 
approach which has been devised in previous articles This approach 

can be seen as a very efficient way to obtain kinetic energy operators for the family of 
polyspherical coordinates. The formalism can be applied whatever the number of atoms 
and whatever the set of vectors: Jacobi, Radau, valence, satellite, etc. The formalism is not 
restricted to total J=0 and hence the operator may include overall rotation and Coriolis 
coupling. 

The protonated water dimer system is described by six Jacobi vectors. The choice of the 
Jacobi vectors is not unique and several clustering schemes are possible, one natural choice 
for HsO^ is given in Figure [H 

FIGURE m AROUND HERE 

The polyspherical approach straightforwardly provides the expression of the kinetic en- 
ergy operator in terms of the angular momenta associated with the vectors describing the 
system. Here we use the technique of "separation into two subsystems" which is described 
(see Eq. (37) there). The full kinetic energy operator reads (we use the notation 



in Ref. 
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dx = and = throughout the paper): 



f _f Kl^m^ (Jl ■ /+ (La + Lb + if - 2{La + Ln + I) ■ J) J,, 

I v^/ ^ f)2 ID \ , {L'a + L\^- LiA- 2La- Lia)bfa , i^iA' Lia)bfa 
+ l^^~o7r~R~^R^AR^'^) + ^7 — ^2 + ^7 — ^2 

, v^/ D N , {L% + L\^-Lib-2Lb-Lib)bfb , i^iB ' ^ib)bfb 

2^^~9n R "RiB^^B) + ^7 ^2 + ^7 ^2 



i=l 



2m r ^ 2mr'^ 



(6) 



with 
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R, RiA, R2yi5 Rib, R2B, r, are the lengths of the vectors R, Ria-, R2A, Rib, R2B, "f^i 
respectively. 

mo2mH ,, VTLH ,, mo+2mH „„J ™ 2mg(mo+2mH) 

LiA{B) is the angular momentum associated with Ria{b)- 

L2A{B) is the angular momentum associated with R2A(b) and La{b) (= -^1 a(b)+-^2A(s)) 
is the total angular momentum associated with monomer A (B). 

/ is the angular momentum associated with the proton. 

J is the total angular momentum of the system. J = I + La + Lb + Lr and Lr is 
the angular momentum associated with R. 

E2 is the frame resulting from the two first Euler rotations starting from the space 
fixed frame, such that the z^"^ axis lies parallel to R. In other words, the two first 
Euler angles a and [3 are identical to the two spherical angles of R in the space fixed 
frame. 

BFA (B) is the body fixed frame associated with monomer A (B). The two first Euler 
angles aA{B) and I3a(b) of monomer A (B) are defined such as z^^^^^^ lies parallel 
to R2A(b) (in other words they are identical to the two spherical angles of R2A(b) in 
the E2 frame). The third Euler angle 7a(b) is defined such as Ria{b) remains parallel 
to the {{xz)^^'^^^\ x^^^^^^ > 0) half plane. This definition of BFA(B) is similar to 



the one recently used to describe the water dimer (see e.g. Ref 28|]). Note however 
that we have changed the definition of the z^^^^^^ axis which is now parallel to the 
vector between the two hydrogen atoms and not to the vector joining the center of 



mass of H2 to the oxygen atom as in Ref. 28|]. This change was performed to avoid 



the singularity which appears in the KEO when R and z^^"^^^^ are parallel. 

It should be emphasized that all the angular momenta appearing in the kinetic energy 
operator are all computed in the space fixed frame but projected onto the axes of different 
frames. This last point, as well as the properties of the corresponding projections, is 



addressed in detail in Ref. 29|]. We recall here only the main point, i.e. the fact that if a 



vector is not involved in the definition of a frame F, the expression of the projections of 
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the corresponding angular momentum onto the F-axes in terms of the coordinates in this 
frame, is identical to the usual one in a space fixed frame. This very useful property will be 
utilized several times in the following for instance to obtain Eqs. ( I13fl4lll5f25f28ll29p . 



We now slightly change the volume element for the lengths of the vectors (as usual), and 
assume J=0. The operator can then be recast in the following form: 



^ 2iiR^' 2/iRi?2 2m ^ 2mr2 

j=l 



. V^/ _^a2 N , i^B + ^iB • - 2Lb ■ Lib) {L[b • -^ig) 

Z^iB Zfl2B^2B ZHiBlt^^B 



BFB 



i=l 



which is to be used with the volume element 



(7) 



dRdRiAdR2AdR\BdR2BdrdaA sin (^Adf^Ad'jAdaB sin f^Bdf^Bd'jB sin OdOdifu sin OiAdOiA sin OiBdOiB 

(8) 

with 

< R, RiA, R2A, Rib, R2B, r < 00 
< f3A,f3B,0,9,A,0iB < vr 

< aA,'yA,OiB,lB,fH < 27r 

We have 16 degrees of freedom instead of 15 since we have not defined the third Euler 
angle yet. The angles are depicted in Figure [2j 

FIGURE M AROUND HERE 

• 6 and (fn are the spherical angles of the proton in the E2 frame: 6 is the angle between 
z^"^ and the vector r, ipn describes the rotation of r around z^"^. 

• Pa and are the spherical angles of R2A in the E2 frame. The angle Pa describes 
the rocking motion of water A fragment while aA describes the rotation of the water 
A fragment around the R = ze2 axis (the same for B). 
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• 9iA and 7^ are the spherical angles of Ria in the E2A frame, the E2A frame is the 
frame obtained after the two Euler rotations ua and 13a with respect to the E2 frame. 
Consequently, 6ia is the angle between Ria and R2A and 7^ describes the rotation 
of RiA around R2a- Hence 7^ describes the wagging motion of the water A fragment 
(the same for B). The fact that we have used the separation into two subsystems 
for the two water monomers (see Eq. (37) in Ref. 15|) explains why the angles of 
RiA{B) are defined with respect to R2A{b) ^i-nd not with respect to R as in the standard 
formulation of the polyspherical approach (see Eq. (65) in Ref. [l^). This allows us 
to have purely intramonomer angles instead of angles mixing the intramonomer and 
intermonomer motions. 

Since {.L\b)) bpa^b) = ^^\b)) e2 ^^^^^ ^^(s) ^^^^^^ ^\{b)^a{b) since 

the projections of La{b) onto the axes of the E2 or BFA (B) frames are hermitian: see Eqs. 
( IT3]ll5|25ll29l) below), we can rewrite the operator as: 



1=1 1=1 

( A^BFAi2^^Jl2 2h2aR2A ^ ^^^^^^2fIjiR^ 2fi2BRlB^ 

+ {L\^ ■ Lia)bfA^17:, ^ + B2~) + (^IB • ^lB)BFBi7r-, ^ + TT^ B2~. 

ZfliAK^A ^^'2AJ^2A 'ifJ'lB^iB '^fJ'2B^2B 

_ {La ■ Lia)bfa _ (Lb ■ Lib)bfb 

fl2AR2A 1^2bR\b 

'7i 



{La ■ Lb)b2 , )e2 _l (-^a + Lb) ■ Ie2 



(9) 



This operator could be straightforwardly used along with an adequate primitive basis 
set of spherical harmonics and Wigner rotation-matrix elements. However, in this work, 
we prefer to employ a direct product primitive basis such as a DVR for all the degrees of 
freedom (only because the latter basis set is numerically more efficient for our problem). We 
have then to explicitly express of the angular momenta in terms of the coordinates, i.e. we 
have to detail the following terms: 

• (1) iLA)BFA^ {Lb)bFB^ (LlA ' Lia)bfA^ {L\b ■ Lib)bfb 
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• (2) {La ■ Lia)bfa^ (^b ■ Lib)bfb 

• (3) {La-Lb)e2 

• (4) 4 

• (5) {La + Lb) ■ Ie2 

but before doing so we have to further specify the coordinates: 

• (i) to avoid the singularity we use the BF Cartesian coordinates x, and z for the 
proton (the definition of BF frame is given in (ii)). Hence we use, d^pSF = xdy — yd^- 

• (ii) we shghtly change the coordinate system (from E2 to BF) and exphcitly define the 
third Euler angle 7 of the system and then the BF frame. We define 7 such that it is 
identical to the first Euler angle of monomer A (here, we slightly break the symmetry 
between the two monomers). Consequently, we have: 



7 = «A 



a = aB — a. A 

,„BF ,„E2 



(10) 

since J=0, we can then put: 



= 



(11) 



and find 



do^A = -9a - d^BF = -da - {xdy - yd^) 

(12) 



Note that this transformation does not affect the other angles. 
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(iii) in order to have hermitian conjugate momenta for the angles, we use U/5^ = cos f^A, 

= cosPb, Vl0^a = COS 9 1 A, ne^g = cos6'ib. 
For each angle 77, we have (with u^ = cos?]): (9^ = — sin rjdu,^ = — ^1 — u'^du^ and 

9j = -a.„,/r 



u 



2 



We have now 15 degrees of freedom: 

< R, RiAi R2A, Rib, R2B < 00 

-1 < Ui3^,upg,ue,A,ue,g < 1 

< a, 7a, 7b < 27r 

—00 < x,y,z < 00 

Let us now detail all the terms appearing in the kinetic energy operator: 

• (1) For Lbbfei the situation is simple since monomer B is not involved in the definition 
of the body fixed frame of the whole system. Thus, we have : 

LbxBfb = ih^^das + ihsiwyB smpBdu^^ - ih^^cos-fBd^g 

LByBFB = -ifi^^daB + ih cos -fB sin PBdufSB + ^^ill!^ sin 7^9^^ (13) 

LbzBfb = — ihdyg 

and 

4 BFB = -f^\d^,B (1 - ^L)^'^.. + T^T^(C + - '^upud.M) (14) 

PB 

Similarly L\ ^pj^ can be seen as the total angular momentum of monomer A projected 
onto the axes of the Body fixed frame of the monomer. We obtain the usual expression: 

LaxBfa = ih^l^da^ + ihsin-fAsin(3Adu0j, - ih^^cos-fAd^^ 

LAyBFA = -ifi^;^daA + ihcos'jAsinpAdu^A + ih^^ sin 'jAd-y^ (15) 



LazBfa = — ihd^y^ 



and 



L'abFA = -^\du,^{l - u'^Jdu^^ + _^ 2 (^L + ^7A - 2M/3^^«a^7a)) (16) 



12 



However, since monomer A is involved in the definition of tlie tfiird Euler angle of the 
system, we apply the change of coordinates, Eqs. ffTOj) and f|T2l) : 

{dl + x^d'y + y'dl - xd^dyv - ydyd^x 
+ 2xdydc, - 2yd^da + <9^^ + 2up^d^d^^ 
+ 2up^d^^xdy - 2up^d^^ydx) 

(17) 

For the projections of Ria{b) onto to the BEA (B) axes, the situation is not straight- 
forward since Ria{b) is involved in the definition of BEA (B). However, since we have 
used the same convention for the BE frames as in our previous articles^the expressions 



of these projections is already known (identical to Eq. (47a) in Ref. l^] for instance): 



LlAyBFA = ihsaidxAdug^A (^^) 

Lij^^BFA = — ihd^^ 

(This expression can be obtained starting from the E2A frame which is the frame 
resulting from the two first Euler rotations ^a and Pa- The expression of the projec- 
tions of LiA onto this E2A frame are the usual ones since Ria is not involved in the 
definition of this frame. Applying then the third Euler rotation and the coordinate 
transformation leading to the BEA coordinates yields Eq. fITS]) ). Note that Li^^sfa 
is not hermitian (this is due to the fact that Ria is involved in the definition of BEA). 
However, at the end we obtain the usual expression: 



{L\a ■ Lia)bfa = -^'(^n,,,(l - uljd^,,, + ^4^5^) (19) 



and exactly the same for B: 



(lIb ■ Lib)sfb = -^'{d^e.si^ - + Y^-^d'J (20) 

1 LLn _ 
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(2) The terms in {La{b) ■ are not obvious since Ria{b) is involved in the 

definition of BFA (B). First, let us explicitly highlight the hermiticity of the term: 



(La--^ia)^^^ = 1/2{D^ ■ LiA + L\ji^ ■ La)bfa 

(21) 

Combining Eq. ( JTSl) and Eq. ( JTSl) (as well as Eq. ( JTOl) and Eq. ( fT2l) ). we obtain: 
{La ■ Lia)bfa = 

+ 29^ + (a„ + a;9j^ - yd^ f^^^f sin 6'ia9„ 

+ sin /^A cos 7a sin 6'iA9„g^^ + d^^ sin 7a sin OiAd^g^^ 

Sin 

sm uiA sm Pa sm uia 

+ "^"si 4 sin OiAida + xdy - yd^,) 

sm Pa 

+ a„ sin /3a cos 7a sin 6'ia9„ + (9„ sin 7a sin OiaO^^] 



(23) 

and for B (the situation is a little bit simpler since monomer B is not involved in the 
definition of the third Euler angle see Eq. f|T2|) ): 

{Lb ■ Lib)bfb = 

^eiB cos 7b ue^^ ■ n . o ^'^is ""fe a 

~ T .in^ .in/3 ^-^^ + ^"''s — ^ ^^^^B - ^d^B COS 7B —75-^73 

z sm f7i5 sm Pb ^smfifj smfifj smp^j 

+ 2d^ - dj^-^ sin 6'ib9„ + 9„ sin /3b cos 7^ sin 6'ib9„ + ."^^ sin 7^ sin 6'ib9„ 
sm /3b '^^ sm /)b 

Mg^^ cos 7b Mgjg • /3 • Q 

smt^iBsm /3b smfc'iB ^ 

sin 6'ib9c, + 9„ sin /3b cos 7b sin 6'ib<9„ 

i^sm/^B 

<9„. ^ sin 7b sin ^ib^^s] 

IS sm /3b 

(24) 



- Que 
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• (3) For the term in 1/2{L\ ■ Lb)^2^ note that the expression of the projections of 
La{b) onto the E2 can be seen as those of total angular momentum of a system onto 
the axes of the Space fixed frame. Indeed, since monomers A and B are not involved 
in the definition of E2, we have (similar to Eq. (18) in Ref. [l^): 



La^e2 = ih cos a A-^^da^ - z/lsina^ sin/3A<9„^^ - ifi^^d^A 
LAyE2 = ihsinaA^ij^da^ + ihcosaA sin pAdu^^ - ^^i^^7A (2^) 
LazE2 = —ihda^ 
After, noticing that 

{La-Lb)e2 = T^i^A' + L^B ' ^a)e2 



we obtain after using Eq. f[T^ : 



(26) 



-[L\ ■ Lb)e2 = - — [-oleosa—— — r^<9a 
2 '^'^ 2 sm/?A sm. 13b 

- {xdy - yd^) cosa-^ — ^-^^^ 

smp^smpe 

+ ^7A • n + da sin a sin/3ij5„ + {xdy - yda^) sin a sin/^s^^^ 

smp^smps smp^ '^^ smp^ '^■^ 

+ 9„ cos — 5~^7i3 + sin a sin /^A ."^^ da 

smpAsmpB ^ sm Pb 

+ {xdy - yd^) cosa-^ — ^—5-^73 
smp^smpe 

+ 9„ cosasin/^Asin/^Bc}^ - ^ug, sin a^^^^^S^^ 

PA PB PA sm Pb 

- dy^ cosa^^- — r^-da + dy^ sina ^^^^^ g^ - dl - da{xdy - yd^,)] 

smp^smps smpA 

(27) 

(4) Since the proton is not involved in the definition of the whole body fixed frame, 
the expression of the projections of / onto the E2 (or BP) axes is the usual one. At 
the end, we obtain: 



(OL = (OL = -^\y'dl + z^dl + z^dl + x^dl + x^dl + y^dl 

-ydyd^z - dyyzdz - d^xzd^ - xd^d^z - xd^dyy - d^xydy) (28) 
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• (5) For the last term {La + Lb) • ~ + ^bV ■ I + 1^ ■ {La + Lb)bf^ need to 
know the expression of the projections of La{b) onto the body fixed axes. For Lb, 
there is no particular problem since B is not involved in the definition of 7 and thus 
of the BF frame: 



LorBF = ih "'^f COS ada — ihsin a sin PBdu .r, — ih 

LrvBf = ihsin a "^f da + ih cos asm dBdu„r, — ih 'j'™" d^„ (29) 

LbzBF = —ihda 

— * 

For La, it is less straightforward. However, applying the third Euler rotation to Eq. 
fl25|) and using Eq. ([12]) yields: 



'AxBF 



■_UiA_ 

'sin /3a 



{da + xdy - yd^) 



sin 13 A A 



LAyBF = i sin (SaOu^^ 
LazBf = i{da + xdy - ydx) 



(30) 



and thus: 

{La + Lb) 



' ^E2 = T^ii^A + Lb 



2 

sin Pb 
9 



^U2^^d -2 ^ d 



2 smp^i smp^i 
sinasinfjBOupB 



— -— cos aOa + oa -— cos a 

mpR sm Pb 



- dupB sin a sin Pb - 2 ^^^^ d^^ ) {yd^ - zdy) 

sm 



0111 L. 

+ {sin pAduA + dupA sin/^A 

+ sma — ac + aaSma — V cos asm pBOuf,,, 

smpB sm Pb 

+ dupB cos a sin Pb - 2^^^^d^g){zdx - xd^) 

sm Pb 

+ 2{x'^dy + y'^dl - xd^dyy - d^xydy)] 

- h^[ .^t, i^dyyd^ + xydyd^ - 2xdyZ - 2d^y'^d^ + d^ydyZ + d^dyyz)] 
Sill Pj\ 

(31) 



The final expression of the operator in Eq. ([9]) is recast as T = Ti + T2 + ^3 + ^4 with 
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2 A2 /'72> 



^£2 



+ ^^^^^^^^2;;!^ + 2;;!^^ + ^^^^^^^ 



(32) 



^ {La-Lia)bfa {Lb-Lib)bfb 

lo — - 



2B 

(33) 



T3 = i^^Al^ (34) 



(35) 



The expression of the operator in terms of the 15 degrees of freedom: 
-R, R\A^ R2A, Rib, R2B, X, y, z, a, u^^,'yA, upg,'~^B, ug-^A^ ^fis reads: (note that the hermiticity 
clearly appears) 



fi = (-- ^4) + y (-- ^4 J + y (-- ^4 )- -^(52 + 92 + 92 

^ 2fiR j^"- 2fi,A j^^ 2fi,B 2m^ ^ y ' 

- f^\du,^ (1 - ^L)^-^4 + 1 \ + + ^'-92 - - ydyd^x 



+ 2M/3^^7Aa;5|/ - 2M/3a^7A^^x))( ^,/p2 + o,. d2 



/^^(^«.. (1 - ^L)^^.. + Y^f. + C - ^MM^2 + 



- ^^(a.,,,(i-<ja„,,, + -^9^j( 



(36) 
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+ 



+ 



Zfi2A^2A sm UiA sm Pa 

'""A ^'"^ ^^^-^-^ " ^^-^^ """^ ^"^^^^^^ 

^ sm (7i^ sm sm Uia 

2dl + (a« + {xdy - sin 6'ia<9„ 

sm 1^ 

+ duf}^ sin /5a cos 7a sin 6'iA<9n9^^ + <97^ sin 7a sin 6'iA9„g^^ 

sm fj 

sm uiA sm Pa 

+ c^7A^-%^sin/3Asin7A<9n„ 
sm OiA ^ 

+ "^^s, . ^7^^ sin 6'ia(9„ + {xdy - ydn,)) 
sm Pa 

+ 9„ sin/5ACOS7Asin6'iA9„ -^^^sin7Asin6'iA9^4] 



~ ^'^7S^-Z COS7ij . Cl^g 

smt^iB smpB 



+ 2(9^ - d j'^^'^f sin 6'ib9„ + (9„ sin /3b cos 7^ sin 6'ib9„ 
sm/3_B 

+ c^7s^^^^sin7Bsin6'iB9„ 
sm/^B 

— — sm [3b sm 7Bd«.^ 

3m U^ R '^■s 



I ^eiB COS 7b ^eiB 

^ sin ^iB sin /5b sin 9^^ 

iin 7b 

— — sin diBda + 9„ sin (3b cos 7b sin 6'ib9„ 
im/^B ^ 



- ^— sm^^iB 

IS sm /5b 

+ <9„ ^4^^sin7Bsin6'iBC?7B 
i^sm/^B 
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z smpAsmpB 

- {xdy - yd^)cosa-^ — r^^^a 

sm Pa sm pb 

- da cos — 'r^i.xdy - yd^) 

sm Pa sm Pb 

+ ^C'ta ^-5 — ^-iT^^B + sm a . sm /^sd^.^ 
smp^smpe smp^ ^ 

+ 9„ sin a sin /?b(9„ + {xdy - yd^) sin a sin PbOu 

smpA smp^ 

+ dup^ sin a^^^ sin ps^xdy - yd^) 

+ (9^ cos «^-5 — + ^^75 cos — — 
sm Pa sm p ^ sm p^ sm p ^ 

+ (9„. sin a sin /3a ^^^^c^Q + 9^ sin a sin /3a c?^ 
'^'^ sm/5B sm/^B 

sin /3a sin /^b^^^^ 



-(9 



2{xdy - yd^) cosa-^ — r^— 
sm/?Asm/?B 

+ duf,^ cos a sin /3a sin hdufj^ + c^n^^ cos a 



. sin/3A„ ^ . sin/^Ajj 

a^. sm a —o^„ — sm a —Ou» 

sin /3b sin /3b "-^ 

COS a^— — - c'q cos tt^— — — 5-a 
sm /3a sm p b sm /3a sm p b 

sin /3d . sin /3d . 



„ . 0111 fJB o ,0 . ^m p, 

+ sma . + d„ sma^— - 

sm /3a ^ ^ sm /3a 

- 2dl-2da{xdy-yd^)] 



d. 



lA 



19 



-ydydzZ - dyyzdz - d^xzdz - xd^d^z - xd^dyy - d^xydy) 



2/1/2 i?2 sin /^A sin (5 a 

+ — cos aOa + Oa TT- COS « — sm « sm pBOutiB ~ ^u,3B sill ct sm pb — 2 T^-o-yr, 

sinpB sinpB sm Pb 



{yd, - zd. 



yj 



+ {sin PAdupA + ^upA sin/3yi 

+ sma^ — --da + da sma^ — — + cosasmpBOu^B + cosasm/Js - 1- — -^o^g) 
sm p B sm Pb sm p £j 

+ 2(x29j + 1/252 - xd^dyy - d^xydy)] 

~ n ^ T^y l {xdyyd, + xydyd, - 2xdlz - 2dr,y^d, + d^ydyZ + d^dyyz)] 
ZjjifiK sm Pa 

(40) 

This completes the derivation of the kinetic energy operator of HsO^. We em- 
phasize again that this operator is exact. Furthermore, the correctness of the deriva- 
tion and implementation of the KEO was checked by comparing it with numerical 
results provided by the program TNUM 30|. A KEO can formally be written as 
T = (1/2) Gij{i^)didj + FM)dj + Ve^tra{q). TNUM computes G, F and Ve^tra nu- 
merically. We have checked that the numerical values of all the functions Gij{q) at several 
non-symmetrical grid points q agree with those provided by the program TNUM. The func- 
tions Fj{q) are determined through the hermiticity of the KEO. As our KEO is obviusly 
hermitian, there is no need to check the Fj. 



B. Hierarchical Representation of the Potential using Mode-Combination 

The exact PES {v{q)) for B.5O2 is a function of the 15 internal coordinates previously 
defined, where "exact" refers to the full dimensional PES of Bowman and collaborators 
a. The calculation of vltaational levels o. the IR spect.u. .cu.es a high accuracy . the 
representation of the Hamiltonian. The exact and trivial representation of v{q) on a product 
grid is unfortunately beyond current computational capabilities: in this case the potential 
would be given on a grid of ~ 10^^ points (~ lO'^ TB of disk space). A direct use of the 
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potential is hence impossible and in particular it is impossible to convert the potential to 



31 



321 ] since the full product grid is 



MCTDH product form by using the potfit algorithm 
needed for such a transformation. The problem of representing a high- dimensional PES for 
quantum-dynamical computation has already been considered in the context of Multimode 
simulations 18|, as well as in the more general context of the high dimensional model 



representation (HDMR) 
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33 



34 



35|. 



In such hierarchical representations a multidimensional function dependent on / variables 
is approximated as: 



/ / / 

a=l a<f3 a<l3<'y 

where q is a coordinate vector and v{q) denotes the hierarchical approximation to v 
component functions in Eq. fj^T]) can be determined by minimizing the functional 



(41) 



he 



35| 



b(q) -v{q)fw{q)dq, 



(42) 



D 



where w{q) is some weight function on the integration domain that determines the form of 
the c omp onent functions. Taking w{q) = Yli'^Kign) leads to component functions of the 



form 



34 



35 



36| 



,(0) 



,(0) 



(43a) 
(43b) 



f 



qp) = [ n wMv{ci) rfq"^ - v^^\q^) - v^^\q^) - (43c) 



There q"" represents a vector of coordinates in which the a-th component has been removed, 
qO/3 represents a vector in which a and (3 components have been removed, and so on. Non- 
separable weights have also been considered, which lead to more complicated expressions for 
component functions (see for example the Appendix in Ref. 1351 ). It is possible to evaluate 



the integrals in Eq. 
called RS-HDMR 33| 



34 



35 



random sampling of the coordinate space, leading to the so 



361]. However, the cumbersome multi-dimensional integrals 
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can be solved trivially by choosing the particular weight 



n ^(^" - 



a=l 



which corresponds to the cut-HDMR approximation 
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(44) 



There a„ is the a- 



th component of point a, the reference expansion point in coordinate space. Using the 

jwhich we will also refer to as 



definition of w{q) in Eq. (jH]), the different terms in Eq 
uncombined clusters (UC), are given up to second order by 
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,(0) 



Via 



,(0) 



af3\ 



,(0) 



(45a) 
(45b) 
(45c) 



where higher orders follow trivially. Cut-HDMR representations of a PES - also called n- 
mode representation - have been already used successfully to accurately compute vibrational 



energy-levels of molecular systems 



37( 1 and reaction rates for molecule-surface scattering 



. Unfortunately, the use of a HDMR representation of the PES of the form of Eq. (HT]1 
leads to a combinatorial increase in the number of terms as the order of correlation increases. 
At the same time, higher order terms are given on grids with an exponentially increasing 
number of points, which leads to quick stagnation in the maximum correlation order that 
can be practically included. 

Instead of directly adopting the expression in Eq. (14T|) for v{q) we make use of the fact 
that the wavefunction in Eq. ([T]) is given in terms of combined modes. Hence we define the 
potential also as a function of the combined modes. The use of combined modes instead 
of coordinates as base of the hierarchical expansion attenuates the combinatorial increase 
in complexity found when using Eq. pTl) while still retaining the simple evaluation of the 
expansion terms given by Eq. (145|) and the inclusion of high-order correlations. 
For / coordinates q = p particles or combined modes Q = [Qi, 



are defined, such that Qi = [qi\ ■ ■ ■ , q^f'] and Yl'i=ifi ~ /• '^^^ reference potential 
be equivalently given as a function of the combined modes, V^(Q) = v{q). The general 
hierarchical expansion of Eq. pTl) is now written in terms of the combined modes Qi, 
instead of coordinates qa, which up to second order reads: 



can 



V{Q) = Vo + Y, V}'\Q^) + ^^f 



(46) 
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First order V^^^ and second order V^p combined clusters (CC) are defined in an analogous 
way to Eq. (145|) : 

= V{s.) (47a) 
V|'\Q^) = ViQf,^')-Vo (47b) 
V!^\Q^, Q,) = V{Q,, Qf, a^^) - v}'\Q,) - vj'\Q,) - Vo (47c) 

First order CC can be given directly in multidimensional product-grids which are direct 
products of the corresponding ID DVR grids since the wavepacket is expanded in sums 
of products of SPFs which are defined on the same combined grids as used in the cluster 
expansion. Second order and higher CC are conveniently and accurately represented as sums 



of products of multidimensional mode- functions by employing the potfit-algorithm [31|, |32 |. 

On the other hand, first order CC can be exactly represented as a cut-HDMR expansion of 
the form of Eq. fHT]) up to order /j in terms of the UC made of the coordinates [qi \ ^2*^ • • •] 
in Qi, while analogously, second order CC can be given by a cut-HDMR expansion up to 
order fi + fj in terms of the corresponding UC. When inspecting the expansion in Eq. (146|) 
one realizes that the correlation between three coordinates belonging to different modes, 
i.e. [qa\qj^\q\''''], is accounted for up to second order with respect to the uncombined 
representation since only the UC v^^^{qa \ q^p), Va]{qa \ q^"^) and v^p^{qP , g^^^) are contained 
in the second-order CC V^p{Qi, Qj), V^f\Qi, Qk) and vj^\Qj, Qk) respectively. Thus, using 
a second order expansion of CC one implicitely introduces a cluster-selection scheme of higher 
(up to fi + fj) order UC. The higher-order cluster-selection scheme is determined by how 
coordinates are grouped together into combined modes. The same reasoning can be extended 
to higher than second order CC. By combining coordinates which are strongly coupled it is 
then possible to get a high-accuracy potential while still using a reduced number of CC. 

FIGURE |3] AROUND HERE 

Assuming that the grid representation of each coordinate uses on average N points, the 
number of points in coordinate-space, Ntoti used to define the UC, i.e. without mode com- 
bination, is given by 

^- = E(«)^"' (48) 

a=0 ^ ^ 

where h is the maximum allowed order for the clusters and / is the number of degrees of 
freedom. In case the coordinates are combined into modes with m coordinates each, the 

23 



number of grid points needed to represent the clusters is given by 

Ntot = Y.\ "^W"^. (49) 



j=0 ^ ^ 

The number of points for different values of the parameters defining different clustering 
schemes is depicted in Fig. [31 The horizontal axis represents /i, the maximum order of 
clustering in terms of the uncombined coordinates, i.e., the maximum order of the UC in 
the expansion. A total of / = 15 coordinates is assumed. A horizontal line is drawn at 
10^, which, tentatively, constitutes a practical limit to the total number of grid points that 
can be used both concerning the generation of the clusters and their subsequent use in the 
dynamical calculations. The example in Fig. [3] assumes that there are A^ = 10 points per 
coordinate. Using 3D modes it is possible to include 6th order UC with around 10^ grid 
points, which is of the order of the PES presented in section UlI CI The inclusion of up to 5th 
order UC without mode-combination would require more than 10*^ points while the inclusion 
of 6th order UC would require around 10^° points. We emphasize again, however, that the 
representations up to order h with and without mode-combination are not equivalent. The 
representation without mode-combination contains all the possible UC of coordinates up to 
order h. In the case of using CC one is implicitly selecting a subset of UC. A definition of 
meaningful mode combinations should however be possible in most cases. Such a definition 
would be based on chemical common-sense, e.g., coordinates belonging to the same chemical 
group or to the same molecule in a cluster are good candidates to be combined. As a final 
remark, we note that there is ongoing effort by other groups 3^, l35|] to use parametrized 
functions instead of grids to describe the clusters. This approach may turn out to be more 
efficient than a grid representation. 



C. Potential Energy Surface for the H5O2 cation 

In order to construct the PES for the HsO^ cation employing the approach described 
above one must start by defining the combined modes that are going to be used. In the 
present case the following five multidimensional modes are selected: Q\ = [z,a,x,y], Q2 = 
[lA, 1b], Q3 = [R, up^, upg], Qi = [RiA, R2A, U0^a] and Q5 = [Rib, R2B, uq^b]- It is convenient 
that coordinates x, y and a are grouped together due to symmetry conserving reasons 
which are exposed below. Modes Q2 and Q3 contain the wagging and rocking coordinates, 
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respectively. Modes Q4 and contain the Jacobi coordinates which represent the internal 
configuration of each water molecule. Coordinates z and R are good candidates to be 
combined together, as will be discussed in Section llVCi They are not combined here since 
this would require the definition of a 6th mode. We have, in fact, started to do some tests 
with a 6 particle mode-combination, and we give some brief remark on this below. The 
definition of the underlying ID grids is provided in Table [H 

TABLE H AROUND HERE 

Following the procedure outlined above one may select a reference point in coordinate 
space and proceed straightforwardly to the computation of the clusters defining the cluster 
expansion. Instead of this, the PES expansion is defined in terms of M = 10 reference points 
and the weight in Eq. P2l) takes the form 

1 

1=1 

The reference points are located on or very close to stationary points in the lowest energy 
regions of the PES. After Eq. (ISO]) the PES expansion is given by 

1 ^ 

vtom = j^Y.^i{ci). (51) 

1=1 

The Vz(Q) terms are given by Eqs. (l46l) and (l47l) . The specific form of Vz(Q) that has been 
used here is given by 

5 45 

V^(Q) = V!^'^ + 5^ V«(g.) + E E + V,%{z,Q,,Q,l (52) 

i=l i=l j=i+l 

where the modes Qi - ■ -Q^ have been defined above. The Vj''"'* term is the energy at the 
reference geometry /. The Vi]'^ terms are the intra-group potentials obtained by keeping the 
coordinates in other groups at the reference geometry /, while the Vjy terms account for the 
group-group correlations. The potential with up to second-order terms gives already a very 
reasonable description of the system. The VJ\,23 term accounts for three-mode correlations 
between the displacement of the central proton, the distance between both water molecules 
and the angular wagging and rocking motions. Note that the primitive grids in each coor- 
dinate are the same irrespective of the reference point used to expand the potential. This 
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means that the average, Eq. (15T1) . can be carried out before the dynamical calculations by 
summing over all the generated grids of the same coordinates for each reference geometry, 
involving no extra cost for the dynamics. The justification for the multi-reference approach 
lies on the nature of HsO^. The protonated water-dimer is a very floppy system featuring 
several equivalent minima and large amplitude motions that traverse low potential energy 
barriers. Thus, the amount of configurational space available to the system at low vibra- 
tional energies is already large and then it is not well covered by a single reference point. The 
property that, for a single reference point, the PES expansion is exact at the reference point 
and hypersurfaces involving the displacement of up to hm modes is lost after averaging over 
several reference geometries. However, the overall mean error is reduced by the averaging. 

The use of several reference points has a further implication which is related to the sym- 
metry properties of the system Hamiltonian. In the case of HsOj' a possible single reference 
geometry to define the PES expansion would be one of the eight equivalent absolute minima 
which belong to the C2 point groups. This choice, however, breaks the total symmetry of the 
Hamiltonian. We have seen that cut-HDMR is exact at the reference point and all hyper- 
surfaces in which up to hm modes have been displaced from the reference point, where hm is 
the expansion order in terms of CC. The description of the rest of C2 points when one is used 
as a reference is thus not equivalent and the whole symmetry is broken. A possible solution 
would be to use as reference one of the two equivalent stationary points, which lie around 
300 cm~^ above the absolute C2 minima. They correspond to a = 90 or 270 degrees, respec- 
tively. Using only one of the points as reference results in a similar breakage of the total 
symmetry of the Hamiltonian. The same happens again if one of the two equivalent 'D2h 
(a = or 180 deg.) stationary points is used, which are even higher in energy. One should 
note that this is a highly symmetrical, multiminima system in which several permutations of 
identical particles are possible through crossing of low energy conformational barriers. The 



vibrational 



levels of HsO^ can be labeled according to the permutation-inversion symmetry 



group QiQ (39I, |40|, which contains the 'D2d point group as a subgroup, but additionally allows 



to permute the H-atoms of either of the two monomers. The two T>2d and eight equivalent 
C2 geometries which are used as reference points are depicted in Fig. HI The difference 
between the structures at the left and right columns is a 180 degrees rotation of a, or equiv- 
alently the permutation of the two hydrogens of one of the water moieties. The way out of 
the symmetry breakage problem is to use all the structures depicted in Fig. H] as reference 
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points of the cluster expansion. The final potential is given then, as discussed above, by the 
average with respect to all the reference points. By using this set of reference points the 
symmetry of the original PES is maintained. Indeed, it will be maintained as long as, for 
an arbitrary selected reference geometry, all the symmetry equivalent points generated by 
the permutations-inversions of the Qiq group are also considered as reference points. 

FIGURE H AROUND HERE 



IV. RESULTS AND DISCUSSION 
A. Quantum Dynamical Calculations 

The kinetic energy and potential energy operators already discussed are used to compute 
the zero point energy (ZPE) of the system and the corresponding ground-state vibrational- 
wavefunction. The algorithm that implements the computation of eigenvalues and eigenfunc- 
tions of the system Hamiltonian within the MCTDH program is called improved relaxation 
and is described elsewhere 4lj. This algorithm is essentially a multiconfiguration self- 



consistent field approach that takes advantage of the MCTDH machinery. All the reported 



42| 



simulations were performed with the Heidelberg MCTDH package of programs 

The comparison between the largest, converged MCTDH calculation and other reported 
results on the same PES is given in Table HTl As a reference we take the given diffusion Monte 
Carlo (DMC) result which has an associated statistical uncertainty of 5 cm^^. A simple 
normal-modes analysis (NMA) with normal modes constructed from the Hessian matrix 
taken at the C2 minimum yields a ZPE of 12 635 cm~^, only 242 cm~^ larger than the DMC 
result. We believe that this surprisingly good result arises from fortuitous error cancellation. 
As will be discussed below, 3 out of the 15 internal degrees of freedom (7^, 7^ and a) are 
not described by a single-well in the lowest energy region of the potential, while the proton- 
transfer motion (z) features a nearly quartic potential which is strongly coupled to the 
water- water stretching (R). The most comprehensive calculations on the vibrational ground 
state based on a wavefunction approach to date are those of Bowman and collaborators ^ 



using the Multimode program [18|. The vibrational configuration interaction (VCI) results, 
both using the single reference (SR) and reaction path (RP) variants are found in Table 
[m These calculations use a normal-mode based Hamiltonian. They incorporate correlation 
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between the different degrees of freedom due to the cluster expansion of the potential [18j and 
the use of a CI wavefunction. The best reported VCI result for the ZPE lies still 104 cm~^ 
above the DMC value. It is worth to mention that before switching to a Hamiltonian based 
on polyspherical coordinates we tried a Hamiltonian expressed in rectilinear-coordinates and 
obtained results similar to those of Bowman and collaborators. 

TABLE M AROUND HERE 

The MCTDH converged result for the ZPE is given in Table [Tll The obtained value for 
the ZPE is 12 376.3 cm-\ 16.7 cm'^ below the DMC value. Table [ml contains ZPE val- 
ues obtained using an increasing number of configurations. According to these results the 
MCTDH reported values are assumed to be fully converged with respect to the number 
of configurations. The deviation from the DMC result must be attributed to the cluster 
expansion of the potential, Eqs. fl51|52l) . 

TABLE ini AROUND HERE 

We give here some technical details regarding the largest MCTDH calculation with 
10 500 000 configurations, which is probably one of the largest MCTDH calculations per- 
formed to date in terms of required computational resources. The calculation was made 
using the recent parallel version of the MCTDH code which is still under development in 
our group. The calculation was run on a 8-processor, shared-memory machine with Intel 
Itanium2 processors, and used a total amount of 13.5 Gb of main memory. 4 Gb were de- 
voted to storage of the Krylov vectors and Krylov times Hamiltonian vectors of the Davidson 
diagonalization procedure. The rest was needed for the representation of the Hamiltonian 
and the mean fields, the vector of coefficients, SPFs and work arrays. Each step of the 
wavefunction relaxation lasted around 13.5 hours of wall-clock time, and consisted of the 
aforementioned Davidson diagonalization of the system Hamiltonian in the current basis of 
SPFs followed by a 2 fs imaginary-time propagation of the SPFs. The whole procedure was 
iterated until self-consistency of the expansion coefficients and SPFs was reached. 

In contrast, the smallest calculation reported in Table UTTl which used 172 800 configu- 
rations, consumed a total amount of memory of around 263 Mb. Each diagonalization of 
the system Hamiltonian plus imaginary time propagation of the SPFs lasted 10 minutes of 
wall clock time using two processors in parallel on the same Intel Itanium2 cluster. The 
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comparison of the two aforementioned simulations illustrates one of the major strengths of 
MCTDH, namely, the usage of variationally optimal coefficients and orbitals leads to an 
early convergence with respect to the size of the multiconfigurational expansion. Relatively 
good results (only 7.4 cm~^ energy difference between the largest and smallest calculation) 
are already obtained using very moderate computational resources. 

B. Properties of the ground state of the system 

The probability density of the ground state wavefunction with respect to some selected 
coordinates and integration over the remaining coordinates is given in Fig. [51 Fig. [5^ 
shows the density along the proton-transfer coordinate z. The probability density is non 
negligible in a range spanning about 1 bohr. Fig. [5]d depicts the density along the a internal 
rotation coordinate. Along this coordinate the system interconverts between two equivalent 
regions of configurational space. The barrier corresponds to planar configurations of the 
whole system and is about 300 cm"^ high depending on the configuration of the rest of 
coordinates. The system can interconvert between both halves even when in the ground 
vibrational state. The dotted curve in Fig. [5)d depicts the density at a 10 times enlarged 
scale and clearly shows a non vanishing density for a = 0, tt. The splitting state arising 
from the barrier along a has also been computed and the splitting energy has been found to 
be 1 cm~^. The small asymmetry observed in the density in Fig. [5b has been analyzed. It 
arises from the fact that the position of the proton is defined relative to one of the two water 
monomers, monomer A, in order to have 15 internal coordinates (see Eq. flTOl) and related 
text). Consequently, a change in a rotates the central proton so that its relative position to 
monomer A remains unaltered. In contrast, in order for the relative position of the central 
proton and monomer B to remain unaltered during a change in a, coordinates x and y must 
change accordingly. We emphasize that the kinetic operator is still exact and not affected 
by this. However, a, x and y coordinates are defined on non-matching grids. The total 
symmetry between monomers A and B with respect to the central proton is conserved in 
the case of a continuous configurational space instead of a discretized one. Thus, symmetry 
is slightly broken due to discretization of the configurational space. Energetically, the effect 
of this symmetry breakage must be well below 1 cm~^ which is the splitting caused by the 
barrier along a, since a larger perturbation would localize the density on one side of the 
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barrier breaking the double well feature. A perfect symmetry after discretization of the 
coordinates would probably be obtained if one switches to cylindrical coordinates {z, p, 0) 
for the proton and uses identical grids for the angles a and 0. 

Fig. shows the probability density along the wagging coordinates. It consists of four 
equivalent maxima, each of them centered at around ±30 degrees from the planar water 
configuration, so that both water molecules are in pyramidal configuration. The probability 
density corresponding to one of the two water molecules in a planar configuration is however 
quite large and indicates a high probability of exchange between equivalent configurations 
of the system in which the water molecules switch between pyramidal geometries. Each 
of these four density maxima corresponds roughly to one of the C2 equivalent minima on 
the PES. A total of 8 equivalent C2 minima are present since the barrier along coordinate 
a divides the configurational space in two equivalent halves. When both monomers are in 
planar configuration the system interconverts between both and configurations by 
rotation along a. 

FIGURE M AROUND HERE 

Only after the introduction of a Hamiltonian based on polyspherical coordinates a sat- 
isfactory description the HsO^ cation was possible. The use of polyspherical coordinates 
allows for the characterization of the system in terms of well defined stretching, bending, 
rocking, wagging and internal rotation motions, each of which corresponds to a single co- 
ordinate. This fact keeps the correlation between coordinates in the PES relatively small, 
and the relation between different coordinates can be usually understood in simple physical 
terms. The representation of the WF given in these coordinates converges then more quickly 
than a WF constructed from rectilinear coordinates and can be more compact. The price 
to pay, however, is a much more complicated expression for the kinetic energy operator. 

C. Quality of the PES expansion 

In order to illustrate the convergence of the mode-combination based cluster-expansion 
PES introduced in Section [III CI the expectation values of the different terms of the potential 
are calculated for the ground vibrational state. These values are given in cm~^ in Table HVl 
The sum of the first order {'^o\V^^\Qi)\'^o) terms is close to 6800 cm~^, half the ZPE, 
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indicating that they carry the major weight in the description of the PES. The second order 
clusters introduce the missing correlation between modes. They have expectation values 
one order of magnitude smaller than the first order terms with one exception, the matrix 
element arising from the V^'^\Qi,Q2) potential. This can be easily understood by noting 
that modes Qi and Q2 contain coordinates z and R, respectively. These two coordinates are 
strongly correlated and indeed they would be good candidates to be put in the same mode 
in an alternative mode-combination scheme. The only third order term that was introduced 
presents a rather marginal contribution to the potential energy of the system. These values 
prove that the PES representation used is of a good quality and rather well converged with 
respect to the reference PES, at least for the energy domain of interest. The square root of 
the expectation value of the potential squared is depicted in the third column. It is a measure 
of the dispersion around the expectation value and should also ideally vanish. The values 
indicate that the PES representation is good, albeit not yet fully converged. Some more 
terms of higher than second order may be added in the future as computational resources 
allow for it. 

TABLE IV] AROUND HERE 

To finish the discussion it is worth mentioning that we have recently started to test a new 
mode-combination. It is based on 6 modes instead of 5. The definition of the modes in terms 
of the coordinates is as follows: Qi = [z,R], Q2 = [x,y,a], Q3 = [7^7 7b], = [""/Sa' ""fe]' 

= [RiA, R2A, uq-^a] Qg — l^iB, R2B1 ''J'9ib\- beginning we started with 5 modes 

in order to keep the number of configurations as small as possible. This has the drawback 
that one of the combined modes, Qi = [z,a,x,y] is very large, thus its SPFs are harder to 
propagate. However, in the reported calculations we use the parallel version of the MCTDH 
code, and the parts which can be parallelized most efficiently are those involved with the 
vector of coefficients. Thus, a 6-modes scheme seems to be more efficient when the parallel 
code is used. In the 6-modes scheme we include all the first, second and some selected third 
order clusters. The obtained values of the ZPE are still preliminary but located in the region 
of -10 cm~^ with respect to DMC, in good accordance to the results of the 5-mode scheme. 
Such results indicate indeed that the reported simulations are robust with respect to the 
mode-combination scheme. 
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V. SUMMARY AND CONCLUSION 



The protonated water-dimer (HsO^) is studied in its full dimensionality (15D) by 
quantum-dynamical wavefunction-methods, using the MCTDH approach. A set of curvilin- 
ear coordinates is used which accounts for the fiuxional, multi-minima nature of the cation. 
An exact expression for the kinetic energy operator of the system is derived using the poly- 
spherical method. A discussion of the different steps involved in the derivation is given. The 
set of polyspherical coordinates introduced allows the description of the different motions of 
the system, namely proton transfer, water-water stretching, OH stretchings, bendings and 
internal rotation by a single coordinate each which also have clear geometrical meanings in 
terms of angles or distances. The PES used in the calculations is that of Huang et al., the 
most accurate PES available for this system to date. The PES must be represented in a 
numerically and computationally adequate way for the quantum-dynamical simulations on 
a discrete multidimensional grid to be feasible. To this end, a variation of the hierarchical 
cut-HDMR method is presented which takes advantage of the mode-combination strategy 
used to represent the wavefunction. Combined modes, instead of the coordinates, are used 
to define the hierarchical expansion. The expansion of the PES is seen to converge quickly 
with the number of clusters in the expansion, which can be attributed to two factors: 

1. A large amount of the correlation is already captured within the modes, which leads to 
quick convergence of the PES with respect to the number of clusters in the expansion. 

2. The set of polyspherical coordinates used allows for the definition of physically mean- 
ingful modes avoiding at the same time artificial correlations (which are due to inad- 
equate coordinates) to appear both in the potential energy and wave function. 

The ZPE of the system is calculated and the results obtained show an excellent agreement 
with previous DMC calculations on the same PES. The reported ZPE with MCTDH lies 
only 16.7 cm~^ below the reported DMC result. This deviation from the "exact" DMC 
result (statistical error 5 cm~^) must be due to the potential representation, as the KEO 
is exact and as it is shown, the reported ZPE is converged with respect to the number of 
configurations in the MCTDH expansion. A fully converged MCTDH calculation needs a 
rather large amount of computational resources, however, it is shown that reasonably good 
results are obtained with a much smaller configurational space due to the optimality of 
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both the coefficients and the SPFs in the MCTDH expansion. The properties of the ground 
vibrational state are analyzed. The central proton delocalizes on a range of about 1 bohr 
along the water-water axis. The fluxional nature of the system is exemplified in the prob- 
ability density along the wagging 7^ and 7^ coordinates. The 2D space spanned by these 
coordinates presents four density maxima which roughly correspond to the C2 minima of the 
system. In going between different higher density regions the system switches between pyra- 
midal conformations of the water molecules. These conformational changes take place along 
low potential-energy barriers. The system is completely delocalized over these low barriers 
leading to a highly symmetric ground-state wave function. The internal rotation coordinate 
a is seen to divide the configurational space in two equivalent halves. The probability den- 
sity for a = and a = it (planar 1^2/1) is non-negligible for the ground vibrational state. 
The tunneling splitting arising from the barrier along a is computed to be 1 cm~^. The 
convergence of the PES expansion is monitored with respect to the expectation values of the 
potential-energy terms which define it, i.e. by inspecting (\l/o|V^j...|\l/o) and (\l/o|^ij...|^o)^^^- 
A good convergence is observed at second order with respect to the combined modes. The 
only third order term present has a rather marginal contribution to the energy with respect 
to l^&o)- Only after switching to curvilinear coordinates the fluxional motion of the highly 
symmetrical HsO^ cation was correctly accounted for. 

The present paper has focused on the definition of an adequate set of coordinates and 
the derivation of the expression of the kinetic energy using the polyspherical method. A 
convenient way to represent the PES for high-dimensional quantum-dynamical simulation 
has been also discussed. The validity of the Hamiltonian setup has been established by com- 
parison to available DMC results in the literature. Moreover, the properties of the ground 
vibrational state have been investigated. Our study provides a picture of the HsO^ system, 
in which the cluster has to be viewed as highly anharmonic, flexible, multi-minima, cou- 
pled system. We show that a converged quantum-dynamical description of such a complex 
molecular system can still be achieved. The companion paper following this one focuses on 
the infrared spectroscopy and dynamics of HsO^. 
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TABLE I: Definition of the one-dimensional grids. denotes the number of grid points and Xi, 
Xf the location of first and last point. The DVRs are defined in Appendix B of Ref. [2j]. 



Coord. 


N 


Xi 


^/ 


DVR 


z 


27 


-1.8 


1.8 


HO 


a 


21 





27r 


exp 


X 


5 


-0.9 


0.9 


HO 


y 


5 


-0.9 


0.9 


HO 


R 


16 


4.2 


6.5 


HO 




7 


-0.5 


0.5 


sin 




7 


-0.5 


0.5 


sin 


lA 


19 


7r-1.8 


vr-Fl.S 


sin 


IB 


19 


-1.8 


1.8 


sin 


RlA 


9 


0.5 


1.8 


HO 


R2A 


9 


2.2 


3.8 


HO 




7 


-0.5 


0.5 


sin 


Rib 


9 


0.5 


1.8 


HO 


R2B 


9 


2.2 


3.8 


HO 




7 


-0.5 


0.5 


sin 
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TABLE II: Comparison of the zero point energy (ZPE) of the HsO^ cation calculated by various 
approaches on the PES by Huang et. al.0]: diffusion Monte-Carlo (DMC), normal mode analysis 
(harmonic), vibrational CI single reference (VCI-SR) and reaction path (VCI-RP) as published in 



]J and MCTDH results. A denotes the difference to the DMC result. The converged MCTDH 



result is obtained with 10 500 000 configurations. Compare with Table Hill 



Method 


ZPE(cm-i 


) A(cm-i) 


DMC 


12 393 





harmonic 


12 635 


242 


VCI-SR 


12 590 


197 


VCI-RP 


12 497 


104 


MCTDH 


12 376.3 


-16.7 
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TABLE III: Comparison of the zero point energy (ZPE) of the HsO^ cation between different 
MCTDH calculations with ascending number of configurations. The A values are given with 
respect to the diffusion Monte Carlo result, 12 393 cm~^ [l^ . 



SPFs per Mode 


N configs. 


ZPE(cm-i) 


A(cm~i) 


(20,20,12,6,6) 


172 800 


12 383.7 


-9.3 


(35,25,15,8,8) 


840 000 


12 378.5 


-14.5 


(40,40,20,8,8) 


2 048 000 


12 377.8 


-15.2 


(60,40,20,8,8) 


3 072 000 


12 376.7 


-16.3 


(70,50,30,10,10) 


10 500 000 


12 376.3 


-16.7 
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TABLE IV: Expectation value of the different terms of the potential expansion (central column) 
and square root of the expectation value of the potential squared (right column). All energies 
in cm^^. The combined modes read: Qi = [z,a,x,y], Q2 = [jaj^b], Qs = [R^up^jUpg], = 
[RiA,R2A,ue^A\ ai^d Q5 = [RiB-,R2B,ue^B\. 



yd) 


(Qi) 


1293.6 


1807.7 


yd) 


(Q2) 


750.6 


966.9 


yd) 


(Q3) 


171.5 


266.9 


yd) 


(Q4) 


2293.2 


3062.8 


yd) 


(Q5) 


2293.1 


3062.8 


y(2) 


{Qi,Q2) 


-526.9 


1037.2 






- / 0.0 




y(2) 


{Q11Q4) 


-27.5 


231.8 


y(2) 


(Ql, Q4) 


-27.4 


231.7 


y(2) 


(Q2,Q3) 


-10.5 


37.6 


y(2) 


{Q2, Qi) 


-24.7 


117.5 


y(2) 


(Q2,Q5) 


-24.7 


117.9 


y(2) 


{Q3, Q4) 


-18.8 


180.9 


y(2) 


{Q3, Q5) 


-18.8 


180.9 


y(2) 


{Q4, Qb) 


1.2 


9.9 


y(3) 


{z,Q2,Q3) 


1.0 


50.4 
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Figure Captions 

Figure [It Jacobi description of the HsO^ system. The vector R connects the two centers 
of mass of the water monomers. The vector r connects the center of mass of the water dimer 
with the central proton. 

Figure [2} Definition of the angles for the HsO^ system. The angles a a and as describe 
the rotation of the water monomers around the vector R, or equivalently around the z-axis 
of the E2- or BF-frame. 

Figure [3J Number of grid points needed for the representation of the clusters of the PES 
expansion, m is the number of coordinates making a mode. 10 grid points per coordinates 
and 15 coordinates are assumed. A horizontal line is drawn at 10®, which tentatively signals 
the maximum practical number of points both regarding their calculation and the use of the 
grids in the dynamical calculations. 

Figure St Geometries of the 10 reference points used in the PES expansion. The view is 
along the 0-H-O axis. Hence only the closest of the two oxygen and the four hydrogens can 
be seen. The difference between the geometries in the left column and each geometry at 
the right column is a rotation of vr along a. Equivalently, the pairs of structures (a,b), (c,j), 
(e,h), (g,f), (i,d) are related by a permutation of hydrogen atoms of one of the monomers. 
The following coordinates are identical for all reference points: R = 4.70 au, x,y,z = 0, 
Ri{A,B) = 1-07 au, R2{A,B) = 2.98 au, Oia{ib) = 0, up^^^^ = 0. Only coordinates a, ja and 
'-fB differ at the 10 reference points. 

Figure [5} For the ground vibrational state probability density along selected coordinates 
and integration over the rest: probability density along the z proton-transfer coordinate (a), 
along the a internal rotation coordinate (b) and on the 2D space spanned by the wagging 
7a and 7^ coordinates (c). The dotted line in (b) corresponds to a 10 times enlarged scale. 
It indicates that the probability density at a = vr is not vanishing. 
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FIG. 1: Vendrell et. al., Journal of Chemical Physics 
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FIG. 2: Vendrell et. al., Journal of Chemical Physics 
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FIG. 5: Vendrell et. al., Journal of Chemical Physics 
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